##------------------------------------
##
##       REPLICATION CODE
##
##    Presidential Directives
##   in a Resistant Bureaucracy
##
##    Journal Of Public Policy
##
##------------------------------------



rm(list = ls())
options(scipen=100)
options(digits=3)
setwd('~/Dropbox/research/book/implementation/jpp/dataverse')


library(arm)
library(stargazer)


#####  Manuscript


#####      - Figure 1



## issuance decision


tiff('figure1.tiff', width = 7, height = 4, units = 'in', res = 600)


par(mfrow = c(1,2),
    mar = c(4,4.5,3.5,1))

x.max <- 3.8
alpha <- .8
k <- 1 #
sigma <- 1
phi <- 1
b.under <- 2 * phi
d <- function(x) { x^2/(2*sigma) }

## identify when removal threat binds
removal.decision <- function(beta, b, d) (  (beta^2 * alpha * (1 - alpha) ) / (1 - beta * alpha) -
                                      (k / (b + d)) )
beta <- .73
b.bar <- 2 * phi * (beta + 1)

## solve for b
removal.threshold <- uniroot(function(y) removal.decision(beta, y, d(y) ), c(0,6))$root


## 1. left panel

plot(0, xlim = c(0, x.max), ylim = c(-2 ,3),
     main = '',
     col = 'white', cex.main = 1,
     xlab = 'Stakes of Directive (b)',
     ylab = 'President\'s Expected Utility',
     xaxt = 'n', yaxt = 'n')
axis(2, 0, las = 1)
b.tilde <- 2.8
axis(1, c(b.bar, b.under, b.tilde), c(expression(bar(b)), expression(underline(b)), expression(tilde(b))))
axis(3, (removal.threshold + 1), cex.axis = .7, '(removal threshold)', pos = 3.12, tick = F)
axis(3, c(removal.threshold), expression(hat(b)), pos = 3.1, tick = F)
axis(3, removal.threshold, '')

epsilon <- .2
polygon(c(b.bar, b.bar, b.under, b.under), c(-4, 4, 4, -4) , col = 'lightgrey',
        border = NA)
polygon(c(5, 5, b.bar, b.bar), c(-4, 4, 4, -4) , col = 'darkgrey',
        border = NA)
axis(3, c(b.tilde - .2, b.tilde + .2), c(expression(b[L]), expression(b[H])), pos = -2.21)
abline(h = 0, lty = 3)
abline(v = removal.threshold, lty = 3)

## 1. no removal rights
lines(c(b.tilde, b.tilde), c(0,-4), lty = 3)
## from 0 to b.under
curve( 2 * alpha * x - 2 * (1-alpha) * d(x) , 0, b.under, add = T, lwd = 2)
## from b.under to b.max ()
curve( alpha * ( beta^2 * 2*x + 2*beta*(1-beta) * (x - d(x)) - (1 - beta)^2 * 2*d(x) )  -
 (1 - alpha) * 2*d(x), b.under, 5, add = T, lty = 2)
## from b.under to removal threshold (for both types)
curve( alpha * ( beta^2 * 2*x + 2*beta*(1-beta) * (x - d(x)) - (1 - beta)^2 * 2*d(x) )  -
 (1 - alpha) * 2*d(x), b.under, removal.threshold, add = T, lty = 1, lwd = 2)

## 2.removal rights

# between removal threshold and b.bar
curve( alpha * ( beta * 2*x + (1 - beta) * (x - d(x))) +
       (1 - alpha) * ( alpha*beta * (x - d(x) - k) - (1 - beta) * (2*d(x) + k) ),
      removal.threshold, b.bar, add = T, lwd = 1)
# to the right of b.bar
curve( alpha * beta^2 * 2*x +
       alpha*beta*(1-beta) * (x-d(x)) +
       (alpha*beta - alpha^2 * beta^2) * (x-d(x)-k) -
       (1 - 2*alpha*beta + alpha^2 * beta^2) * (2*d(x) + k),
      b.bar, 5, add = T, lwd = 1)



## 2. right panel (legend)


par(mar = c(4,2,2,1))

plot(0, xlim = c(0, x.max), ylim = c(-2 ,3),
     col = 'white',
     bty = 'n', xlab = '', ylab = '',
     xaxt = 'n', yaxt = 'n')
legend('topleft',
       lty = c(1, 2, 1, NA, NA, NA, NA, NA, NA),
       lwd = c(1, 1, 2, NA, NA, NA, NA, NA, NA),
       fill = c(NA, NA, NA, 'white', 'lightgrey', 'darkgrey', 'white', 'white'),
       border = c(NA, NA, NA, 'black', 'black', 'black', 'white', 'white'),
       legend = c('Removal rights', 'No removal rights', 'Both (no difference)',
                  'Self-enforcing', 'Enforceable', 'Unenforceable', '',
                  expression(paste(alpha, '=.8, ',
                                   beta, '= .7, ')),
                  expression(paste(sigma, '=1, ', phi, '=1, ', 'k=1'))),
       cex = 1)


dev.off()




#####      - Figure 2



tiff('figure2.tiff', width = 7, height = 4, units = 'in', res = 600)


par(mfrow = c(1,3),
    mar = c(4,3.5,2,1))

beta.l <- .65
beta.h <- .85
alpha <- .9
alpha.l <- alpha * beta.l
alpha.h <- alpha * beta.h

phi <- 1
k <- 1
sigma <- 1
b.min <-  1 / phi

gap <- 0.035

## 1. left panel

plot(1, xlim = c(.5 + gap, 1 - gap), ylim = c(gap , 1 - gap),
     main = 'No Removal Rights',
     col = 'white', cex.main = 1,
     xlab = expression(paste('Agency Loyalty (', beta, ')')),
     ylab = '',
     xaxt = 'n', yaxt = 'n')
axis(1, c(.52, beta.l, beta.h, 1),
     c(.5, expression(hat(beta)(b[H])), expression(hat(beta)(b[L])), 1))
axis(2, c(0, .5, alpha.l, alpha.h, alpha, 1),
     c(0, .5, expression(alpha * hat(beta)(b[H])), expression(alpha * hat(beta)(b[L])), expression(alpha), 1), las = 1,
     mgp = c(10,.5,0))

curve( x * alpha , 0, 1, add = T, lwd = 3)


abline(v = beta.l, lty = 3, col = 'darkgrey')
abline(v = beta.h, lty = 3, col = 'darkgrey')
abline(h = alpha.l, lty = 3, col = 'darkgrey')
abline(h = alpha.h, lty = 3, col = 'darkgrey')


## 2. middle panel

plot(1, xlim = c(.5 + gap, 1 - gap), ylim = c(gap , 1 - gap),
     main = 'Removal Rights',
     col = 'white', cex.main = 1,
     xlab = expression(paste('Agency Loyalty (', beta, ')')),
     ylab = '', # Compliance Rate
     xaxt = 'n', yaxt = 'n')
axis(1, c(.52, beta.l, beta.h, 1),
     c(.5, expression(hat(beta)(b[H])), expression(hat(beta)(b[L])), 1))
axis(2, c(0, .5, alpha.l, alpha.h, alpha, 1),
     c(0, .5, expression(alpha * hat(beta)(b[H])),
       expression(alpha * hat(beta)(b[L])),
       expression(alpha), 1), las = 1,
     mgp = c(10,.5,0))

polygon(c(beta.l, beta.l, beta.h, beta.h), c(0, 1, 1, 0) , col = 'lightgrey',
        border = NA)
polygon(c(beta.h, beta.h, 1, 1), c(0, 1, 1, 0) , col = 'darkgrey',
        border = NA)

abline(v = beta.l, lty = 3, col = 'darkgrey')
abline(v = beta.h, lty = 3, col = 'darkgrey')
abline(h = alpha.l, lty = 3, col = 'darkgrey')
abline(h = alpha.h, lty = 3, col = 'darkgrey')

curve( x * alpha , 0, beta.l, add = T, lwd = 3)
curve( x * alpha , beta.l, beta.h, add = T, lwd = 1, lty = 2)
lines(c(beta.l, beta.h), c(alpha, alpha))
lines(c(beta.h, 1), c(alpha, alpha), lwd = 3)


## 2. right panel (legend)

par(mar = c(4,.5,2,.5))

plot(0, xlim = c(0, 1), ylim = c(-2 ,3),
     col = 'white',
     bty = 'n', xlab = '', ylab = '',
     xaxt = 'n', yaxt = 'n')

legend('topleft', lty = c(1, 3, 1, NA, NA, NA, NA, NA),
       title = 'Hypothetical directives:',
       legend = c(expression(paste('High-stakes directive (', b[H], ')')),
                  expression(paste('Low-stakes directive (', b[L], ')')),
                  expression(paste('Both (', b[L], ' and ', b[H], ')')),
                  ## 'Removal threat not credible',
                  expression(paste('Removal thrshld. surpassed (', b[H], ' only)')),
                  expression(paste('Removal thrshld. surpassed (', b[L], ' & ', b[H], ')')),
                  '',
                  expression(paste(sigma, '=1, ', phi, '=1, ', 'k=1, ', alpha, '=.9, ')),
                  expression(paste(hat(beta)(b[L]), '=.65, ', hat(beta)(b[H]), '=.85'))),
       lwd = c(1, 1, 3, NA, NA, NA, NA),
       fill = c(NA, NA, NA, 'lightgrey', 'darkgrey', NA, NA, NA),
       border = c(NA, NA, NA, 'black', 'black', NA, NA, NA),
       cex = .85)

dev.off()






#####      - Table 1


## 2x2 table created using LaTeX



#####      - Table 2



x <- read.csv(file = 'replication_data.csv')

## summary stats
z <- subset(x, !is.na(post.before.nprm))
nrow(z) # 11,887
length(unique(z$agency.id)) ## 118
## -------------


## 1. reporting directive

m1 <- glmer(post.before.nprm ~ removal.right +
               (1|year.init) + (1|agency.party) +
              (1|agency.id) +
               agency.loyalty +
             reg.stakes,
            data = z, family = binomial(link = 'logit'),
            control=glmerControl(optimizer="bobyqa",
                                 optCtrl=list(maxfun=2e8)))


m2 <- glmer(post.before.nprm ~ removal.right +
             (1|year.init) + (1|agency.party) +
              (1|agency.id) +
              agency.loyalty +
             agency.loyalty:removal.right +
             reg.stakes,
             data = z, family = binomial(link = 'logit'),
            control=glmerControl(optimizer="bobyqa",
                                 optCtrl=list(maxfun=2e8)))


m3 <- glmer(post.before.nprm ~ removal.right +
             (1|year.init) + (1|agency.party) +
                (1|agency.id) +
               agency.loyalty +
             reg.stakes +
             reg.stakes:removal.right,
         data = z, family = binomial(link = 'logit'),
         control=glmerControl(optimizer="bobyqa",
                                 optCtrl=list(maxfun=2e8)))



## 2. review directive

## subset reviewed regulations
y <- subset(x, x$review == 1)
y <- subset(y, oira.stage.proposed == 1 & stage == 'Proposed Rule')


m4 <- glmer(not.returned ~ agency.loyalty + (1|agency.party) +
               (1|agency.id) +
                   (1|year.init) +
               reg.stakes,
             data = y, family = binomial(link = 'logit'),
            control=glmerControl(optimizer="bobyqa", # Nelder_Mead
                                 optCtrl=list(maxfun=2e8)))


n <- c('Removal Rights',
       'Agency Loyalty',
       'Stakes of Regulation (the index)',
       'Removal Rights $\\times$ Agency Loyalty',
       'Removal Rights $\\times$ Stakes of Regulation')

agency.fe <- c('Agency Intercepts', rep('\\checkmark', 4))
year.fe <- c('Year Intercepts', rep('\\checkmark', 4))
admin.fe <- c('Agency-Party Intercepts', rep('\\checkmark', 4))



stargazer(m1, m2, m3, m4,
          ## type = 'text',
          omit = 'Constant|agency.id|year',
          label = 'tbl:comp',
          title = 'Compliance Models (Logit Coefficients)',
          star.cutoffs = c(.1, .05, .01),
          star.char = c('^\\cdot', '*', '**'),
          covariate.labels = n,
          model.names = F,
          dep.var.labels = rep(c(''), 6),
          add.lines = list(agency.fe, year.fe, admin.fe),
          no.space = T,
          font.size = 'small',
          omit.stat=c("f", 'ser', 'adj.rsq', 'aic', 'bic', 'theta', 'll'))



##        - marginal effects


mean(x$reg.stakes)
x.delta <- sd(x$reg.stakes)
x.min <- min(x$reg.stakes)
x.max <- max(x$reg.stakes)


## removal rights agency
c.0 <- invlogit(fixef(m3)['(Intercept)'] + fixef(m3)['removal.right'] +
                fixef(m3)['reg.stakes'] * x.min +
                fixef(m3)['removal.right:reg.stakes'] * x.min)
c.1 <- invlogit(fixef(m3)['(Intercept)'] + fixef(m3)['removal.right'] +
                fixef(m3)['reg.stakes'] * x.max +
                fixef(m3)['removal.right:reg.stakes'] * x.max)
c.1 - c.0

## no removal rights
c.0 <- invlogit(fixef(m3)['(Intercept)'] +
                fixef(m3)['reg.stakes'] * x.min)
c.1 <- invlogit(fixef(m3)['(Intercept)'] +
                fixef(m3)['reg.stakes'] * x.max)
c.1 - c.0





#####      - Figure 3



tiff('figure3.tiff', width = 7, height = 4, units = 'in', res = 600)


par(mfrow = c(1,3), mar = c(2, 5, 2, 1))


# run RE models above (m2 and m4)


## 1. REPORTING DIRECTIVE

m.report.re <- m2

a <- x[!duplicated(cbind(x$agency.abb, x$party)),
       c('agency.abb', 'agency.loyalty', 'party', 'removal.right')]
b <- data.frame(a.int = ranef(m.report.re)$agency.party[,1],
                agency.party = rownames(ranef(m.report.re)$agency.party),
                se = se.ranef(m.report.re)$agency.party[,1])
b$party <- gsub('.*-(.*)$', '\\1', b$agency.party)
b$agency.abb <- gsub('^(.*)-.*', '\\1', b$agency.party)
a <- merge(a, b, by = c('agency.abb', 'party'))

d <- subset(a, removal.right == 0)
a <- subset(a, removal.right == 1)
alpha.0 <- fixef(m.report.re)['(Intercept)']
weight <- 1.65


alpha.a <- a$a.int
sub.rule2 <- fixef(m.report.re)['reg.stakes'] * .434
error <- a$se
a$prob <- invlogit(alpha.0 + alpha.a + sub.rule2 +
                   fixef(m.report.re)['agency.loyalty'] * a$agency.loyalty +
                   fixef(m.report.re)['removal.right:agency.loyalty'] * a$agency.loyalty +
                   fixef(m.report.re)['removal.right'])
a$lower <- invlogit(alpha.0 + alpha.a + sub.rule2 - error +
                    fixef(m.report.re)['agency.loyalty'] * a$agency.loyalty +
                    fixef(m.report.re)['removal.right:agency.loyalty'] * a$agency.loyalty +
                    fixef(m.report.re)['removal.right'])
a$upper <- invlogit(alpha.0 + alpha.a + sub.rule2 + error +
                    fixef(m.report.re)['agency.loyalty'] * a$agency.loyalty +
                    fixef(m.report.re)['removal.right:agency.loyalty'] * a$agency.loyalty +
                    fixef(m.report.re)['removal.right'])

alpha.a <- d$a.int
error <- d$se
d$prob <- invlogit(alpha.0 + alpha.a + sub.rule2 +
                   fixef(m.report.re)['agency.loyalty'] * d$agency.loyalty)
d$lower <- invlogit(alpha.0 + alpha.a + sub.rule2 - error +
                    fixef(m.report.re)['agency.loyalty'] * d$agency.loyalty)
d$upper <- invlogit(alpha.0 + alpha.a + sub.rule2 + error +
                    fixef(m.report.re)['agency.loyalty'] * d$agency.loyalty)

plot(a$agency.loyalty, a$prob,
     ylim = c(0,1),
     xaxt = 'n',
     xlab = '', cex.lab = 1.25,
     cex = 1.2,
     ylab = 'Compliance Rate',
     main = 'Reporting Directive',
     cex.main = 1.25, font.main = 1)
title(xlab = 'Agency Loyalty', line = .8, cex.lab = 1.25)
curve( invlogit( fixef(m.report.re)['(Intercept)'] + sub.rule2 +
                 (fixef(m.report.re)['agency.loyalty'] + fixef(m.report.re)['removal.right:agency.loyalty']) * x +
                 fixef(m.report.re)['removal.right']) , add = TRUE)
segments(a$agency.loyalty, a$prob, a$agency.loyalty, a$upper)
segments(a$agency.loyalty, a$prob, a$agency.loyalty, a$lower)

points(d$agency.loyalty, d$prob, pch = 20, cex = 1.2)
curve( invlogit( fixef(m.report.re)['(Intercept)'] + sub.rule2 +
                 fixef(m.report.re)['agency.loyalty']*x) , add = TRUE)
segments(d$agency.loyalty, d$prob, d$agency.loyalty, d$upper)
segments(d$agency.loyalty, d$prob, d$agency.loyalty, d$lower)




## 2. REVIEW DIRECTIVE


m.oira.re <- m4

## prepare plot
v <- y[!duplicated(cbind(y$agency.abb, y$party)), c('agency.abb', 'agency.loyalty', 'party')]

r <- data.frame(agency.int = ranef(m.oira.re)$agency.party[,1],
                agency.party = rownames(ranef(m.oira.re)$agency.party),
                agency.se = se.ranef(m.oira.re)$agency.party[,1])
b <- data.frame(a.int = ranef(m.oira.re)$agency.party[,1],
                agency.party = rownames(ranef(m.oira.re)$agency.party),
                se = se.ranef(m.oira.re)$agency.party[,1])
b$party <- gsub('.*-(.*)$', '\\1', b$agency.party)
b$agency.abb <- gsub('^(.*)-.*', '\\1', b$agency.party)
b <- merge(b, r, by = 'agency.party')

w <- merge(v, b, by = c('agency.abb', 'party'))

alpha.0 <- fixef(m.oira.re)['(Intercept)']

weight <- 1.65


sub.rule <- fixef(m.oira.re)['reg.stakes'] * .434
alpha.a <- w$a.int + w$agency.int
error <- w$agency.se
w$prob <- invlogit(alpha.0 + alpha.a + sub.rule +
                   fixef(m.oira.re)['agency.loyalty'] * w$agency.loyalty)
w$lower <- invlogit(alpha.0 + alpha.a + sub.rule - error +
                    fixef(m.oira.re)['agency.loyalty'] * w$agency.loyalty)
w$upper <- invlogit(alpha.0 + alpha.a + sub.rule + error +
                    fixef(m.oira.re)['agency.loyalty'] * w$agency.loyalty)

plot(w$agency.loyalty, w$prob,
     ylim = c(0,1),
     cex = 1.2,
     xaxt = 'n',
     xlab = '', cex.lab = 1.25,
     ylab = 'Compliance Rate',
     main = 'Review Directive',
     cex.main = 1.25, font.main = 1)
title(xlab = 'Agency Loyalty', line = .8, cex.lab = 1.25)
curve( invlogit( fixef(m.oira.re)['(Intercept)'] + sub.rule +
                 fixef(m.oira.re)['agency.loyalty']*x ) , add = TRUE)
segments(w$agency.loyalty, w$prob, w$agency.loyalty, w$upper)
segments(w$agency.loyalty, w$prob, w$agency.loyalty, w$lower)


## 3. legend plot

par(mar = c(1, 1, 2, 1))

plot(w$agency.loyalty, w$prob,
     col = 'white',
     bty = 'n', xlab = '', ylab = '',
     xaxt = 'n', yaxt = 'n')
legend('topleft', legend = c('Removal rights',
                             'No removal rights',
                             '',
                             'Stakes of regulation held',
                             'constant at mean value'),
       pch = c(1, 20, NA, NA, NA),
       cex = 1.2)




dev.off()




#####  Appendix


#####      - Table 3 (summary stats)


options(digits=3)

x <- read.csv(file = 'replication_data.csv')

y <- summary(x[, c('post.before.nprm',
                   'removal.right',
                   'sig',
                   'year.init', 'review',
                   'rfa', 'gov.levels', 'legal.dline', 'unfun.mandate',
                   'agency.loyalty', 'reg.stakes')], na.rm = T)

y <- y[-c(2,5,7),]
y[2,] <- y[4,]
y <- y[-4,]
y <- apply(y, 1, function(x) gsub('Min\\.|Max\\.|Mean|:', '', x))
y <- apply(y, 1, function(x) gsub(' +', '', x))
y <- apply(y, 1, function(x) substr(x, 1, 4))

# correct the mean value here by hand
mean(x$post.before.nprm, na.rm = T) # .44

n <- c('First announced before NPRM',
       'Removal Rights', 'Significant Proposal',
       'Year',
       'Proposal reviewed by OIRA', 'Regulatory Flexibility Analysis Required',
       'Lower Levels of Government Affected', 'Legal Deadline', 'Unfunded Mandate',
       'Agency Loyalty', 'Regulation Stakes')

rownames(y) <- n
std.dev <- sapply(x[, c('post.before.nprm',
                   'removal.right', 'sig',
                   'year.init', 'review',
                   'rfa', 'gov.levels', 'legal.dline', 'unfun.mandate',
                   'agency.loyalty', 'reg.stakes')], sd, na.rm = T)

y <- cbind(y, round(std.dev, 2))
xtable(y)



#####      - Table 4 (list of agencies)


x <- read.csv(file = 'replication_data.csv')

length(unique(x$agency.id)) # 118

y <- x[!duplicated(x$agency.id), c('agency.abb', 'agency.id', 'removal.right')]
z <- read.csv('agencycodes.csv')
z$agency.id <- paste(z$agency.abb, z$bureau.abb, sep = '-')
## corrections
z$agency[z$agency.id == 'DOT-STB'] <- as.character(z$bureau[z$agency.id == 'DOT-STB'])
z$agency[z$agency.id == 'DOE-FERC'] <- as.character(z$bureau[z$agency.id == 'DOE-FERC'])
z$agency[z$agency.id == 'HUD-FHFA'] <- as.character(z$bureau[z$agency.id == 'HUD-FHFA'])
z$agency.id[z$agency.id == 'DOT-STB'] <- 'STB-STB'
z$agency.id[z$agency.id == 'DOE-FERC'] <- 'FERC-FERC'
z$agency.id[z$agency.id == 'HUD-FHFA'] <- 'FHFA-FHFA'

y <- merge(y, z, by = 'agency.id', all.x = T)
y <- y[!duplicated(y$agency.id),]
y <- subset(y, select = c(agency, bureau, removal.right))
y$removal.right <- ifelse(y$removal.right == 1, 'checkmark', '')
y$agency[as.character(y$agency) == as.character(y$bureau)] <- ''
y <- y[order(y$agency, y$bureau),]
rownames(y) <- 1:nrow(y)
xtable(y)






#####      - Table 5 (matched regressions)


library(MatchIt)
library(rgenoud)

x <- read.csv(file = 'replication_data.csv')


sum(is.na(x$unfun.mandate))

table(x$removal.right)

## --- set up genetic matching ---
## y <- subset(x, !is.na(post.before.nprm))
## y <- subset(y, select = c(removal.right, sig, agency.loyalty, rfa,
##                           gov.levels, legal.dline, unfun.mandate,
##                           post.before.nprm, year.init, agency.party,
##                           agency.id, reg.stakes))
## match1 <- matchit(removal.right ~
##                   sig + agency.loyalty + rfa + gov.levels + legal.dline + unfun.mandate,
##                   data = y, method = 'genetic')
## save(match1, file = 'matched_data.RData')
## --------------------------------

load('matched_data.RData')
d1 <- match.data(match1)


m1 <- glmer(post.before.nprm ~ removal.right +
               (1|year.init) + (1|agency.party) +
              (1|agency.id) +
               agency.loyalty +
             reg.stakes,
            data = x, family = binomial(link = 'logit'),
            control=glmerControl(optimizer="bobyqa",
                                 optCtrl=list(maxfun=2e8)))


m2 <- glmer(post.before.nprm ~ removal.right +
             (1|year.init) + (1|agency.party) +
              (1|agency.id) +
              agency.loyalty +
             agency.loyalty:removal.right +
             reg.stakes,
             data = x, family = binomial(link = 'logit'),
            control=glmerControl(optimizer="bobyqa",
                                 optCtrl=list(maxfun=2e8)))


m3 <- glmer(post.before.nprm ~ removal.right +
             (1|year.init) + (1|agency.party) +
                (1|agency.id) +
               agency.loyalty +
             reg.stakes +
             reg.stakes:removal.right,
         data = x, family = binomial(link = 'logit'),
         control=glmerControl(optimizer="bobyqa",
                                 optCtrl=list(maxfun=2e8)))



## update with matching weights
m1.matched <- update(m1, data = d1, weights = d1$weights)
m2.matched <- update(m2, data = d1, weights = d1$weights)
m3.matched <- update(m3, data = d1, weights = d1$weights)


n <- c('Removal Rights', 'Agency Loyalty',
       'Stakes of Regulation (the index)',
       'Removal Rights $\\times$ Agency Loyalty',
       'Removal Rights $\\times$ Stakes of Regulation')
agency.fe <- c('Agency Intercepts', rep('\\checkmark', 3))
year.fe <- c('Year Intercepts', rep('\\checkmark', 3))
admin.fe <- c('Agency-Party Intercepts', rep('\\checkmark', 3))



stargazer(m1.matched, m2.matched, m3.matched,
          ## type = 'text',
          omit = 'Constant|agency.id|year',
          label = 'matched',
          title = 'Compliance Models (MLMs, Matched Data)',
          star.cutoffs = c(.1, .05, .01),
          star.char = c('^\\cdot', '*', '**'),
          covariate.labels = n,
          model.names = F,
          dep.var.labels = rep(c(''), 6),
          add.lines = list(agency.fe, year.fe, admin.fe),
          no.space = T,
          font.size = 'small',
          omit.stat=c("f", 'ser', 'adj.rsq', 'aic', 'bic', 'theta', 'll'))







#####      - Table 6 (regressions w/ additional controls)





x <- read.csv(file = 'replication_data.csv')


m1 <- glmer(post.before.nprm ~ removal.right +
               (1|year.init) + (1|agency.party) +
               pres.senate.seats + net.approval +
              (1|agency.id) +
               agency.loyalty + selin.d1 +
             reg.stakes,
            data = x, family = binomial(link = 'logit'),
            control=glmerControl(optimizer="bobyqa",
                                 optCtrl=list(maxfun=2e8)))


m2 <- glmer(post.before.nprm ~ removal.right +
             (1|year.init) + (1|agency.party) +
              (1|agency.id) +
               pres.senate.seats + net.approval +
              agency.loyalty + selin.d1 +
             agency.loyalty:removal.right +
             reg.stakes,
             data = x, family = binomial(link = 'logit'),
            control=glmerControl(optimizer="bobyqa",
                                 optCtrl=list(maxfun=2e8)))


m3 <- glmer(post.before.nprm ~ removal.right +
             (1|year.init) + (1|agency.party) +
                (1|agency.id) +
               pres.senate.seats + net.approval +
               agency.loyalty + selin.d1 +
             reg.stakes +
             reg.stakes:removal.right,
         data = x, family = binomial(link = 'logit'),
         control=glmerControl(optimizer="bobyqa",
                                 optCtrl=list(maxfun=2e8)))


## 2) review directive

y <- subset(x, review == 1)
y <- subset(y, oira.stage.proposed == 1 & stage == 'Proposed Rule')

m4 <- glmer(not.returned ~ agency.loyalty + (1|agency.party) +
               (1|agency.id) +
                   (1|year.init) +
               pres.senate.seats + net.approval + selin.d1 +
               reg.stakes,
             data = y, family = binomial(link = 'logit'),
            control=glmerControl(optimizer="bobyqa", # Nelder_Mead
                                 optCtrl=list(maxfun=2e8)))


n <- c('Removal Rights',
       'Senate Seats (Pres. Party)',
       'Net Presidential Approval',
       'Agency Loyalty',
       'Agency Independence (Selin D1)',
       'Stakes of Regulation (the index)',
       'Removal Rights $\\times$ Agency Loyalty',
       'Removal Rights $\\times$ Stakes of Regulation')

agency.fe <- c('Agency Intercepts', rep('\\checkmark', 4))
year.fe <- c('Year Intercepts', rep('\\checkmark', 4))
admin.fe <- c('Agency-Party Intercepts', rep('\\checkmark', 4))


stargazer(m1, m2, m3, m4,
          ## type = 'text',
          omit = 'Constant|agency.id|year',
          label = 'tbl:comp',
          title = 'Compliance Models (Logit Coefficients)',
          star.cutoffs = c(.1, .05, .01),
          star.char = c('^\\cdot', '*', '**'),
          covariate.labels = n,
          model.names = F,
          dep.var.labels = rep(c(''), 6),
          add.lines = list(agency.fe, year.fe, admin.fe),
          no.space = T,
          font.size = 'small',
          omit.stat=c("f", 'ser', 'adj.rsq', 'aic', 'bic', 'theta', 'll'))



#####      - Table 7 (agency loyalty checks)



x <- read.csv(file = 'replication_data.csv')

p <- seq(0, .35, .05)
f <- matrix(ncol = 4, nrow = length(p))

for(i in 1:length(p)) {

    agencyL <- quantile(unique(x$ideology1), c(p[i], 1 - p[i]) )[1]
    agencyR <- quantile(unique(x$ideology1), c(p[i], 1 - p[i]) )[2]


    x$ideology.gap[x$party == 'rep'] <- abs(agencyR -
                                                x$ideology1[x$party == 'rep'])
    x$ideology.gap[x$party == 'dem' | x$party == 'dem'] <-
        abs(agencyL - x$ideology1[x$party == 'dem' | x$party == 'dem'])


    x$agency.loyalty.new <- rep(NA, nrow(x))
    x$agency.loyalty.new[x$party == 'rep'] <- abs(max(x$ideology.gap[x$party == 'rep'],
                                                     na.rm = T) -
                                                 x$ideology.gap[x$party == 'rep'])

    x$agency.loyalty.new[x$party == 'dem' | x$party == 'dem'] <-
        abs(max(x$ideology.gap[x$party == 'dem' | x$party == 'dem'], na.rm = T) -
            x$ideology.gap[x$party == 'dem' | x$party == 'dem'])


    m1 <- glmer(post.before.nprm ~ removal.right +
             (1|year.init) + (1|agency.party) +
              (1|agency.id) +
               agency.loyalty.new +
             agency.loyalty.new:removal.right +
             reg.stakes,
             data = x, family = binomial(link = 'logit'),
             control=glmerControl(optimizer="bobyqa",
                                 optCtrl=list(maxfun=2e8)))


    beta <- fixef(m1)['agency.loyalty.new']
    beta.se <- se.fixef(m1)['agency.loyalty.new']
    beta2 <- fixef(m1)['removal.right:agency.loyalty.new']
    beta2.se <- se.fixef(m1)['removal.right:agency.loyalty.new']
    f[i,] <- c(beta, beta.se, beta2, beta2.se)
}

g <- data.frame(dem = p, rep = 1 - p, beta = f[,1], se = f[,2], beta2 = f[,3], beta2.se = f[,4])

xtable(g)





#####      - Figure 4





k <- 1
sigma <- 1
phi <- 1
b.under <- 2 / phi
d <- function(x) { x^2/(2*sigma) }

beta <- .7
b.bar <- 2 * (beta + 1) / phi

a <- seq(.5, .99, .002)
b <- seq(.5, .99, .002)
m <- c()

v <- matrix(nrow = length(b), ncol = length(a))

for(i in 1:length(a)) {
    alpha <- a[i]
    for(j in 1:length(b)) {
        beta <- b[j]
        b.bar <- 2 * (beta + 1) / phi
        removal.decision <- function(beta, b, d) (  (beta^2 * alpha * (1 - alpha) ) /
                                                (1 - beta * alpha) - (k / (b + d)) )
        removal.threshold <- uniroot(function(y) removal.decision(beta, y, d(y) ),
                                     c(0,1000))$root
        m[j] <- ifelse(removal.threshold > b.bar, 3, ifelse(removal.threshold < b.under, 1, 2))
    }
    v[,i] <- m
}

betas <- c()
alphas <- c()
group <- c()
k <- 1

for(j in 1:nrow(v)) {
    m <- v[,j]
    for(i in 1:(length(m) - 1)) {
        if(m[i] != m[i+1]) {
            betas[k] <- b[i]
            alphas[k] <- a[j]
            group[k] <- m[i]
            k <- k+1
        }
    }
}

x <- data.frame(alphas, betas, group)



pdf('scope_conditions.pdf', height = 3.5)

par(mfrow = c(1,2),
    mar = c(3,3,2,1))


plot(0,0,
     xaxt = 'n',
     yaxt = 'n',
     xlab = '',
     ylab = '',
     xlim = c(.55, .95),
     ylim = c(.55, .95))
# x axis
axis(1, at = c(.54, .96), label = c('',''))
axis(1, at = c(.54, .96), label = c(.5, 1), cex.axis = .7, pos = .54, tick = F)
axis(1, label = expression(paste('Pr(Appointee is Loyal), ', alpha)), at = .75,
     tick = F, pos = .52)
# y axis
axis(2, at = c(.54, .96), label = c('',''))
axis(2, at = c(.54, .96), label = c(.5, 1), cex.axis = .7, pos = .54, tick = F)
axis(2, label = expression(paste('Pr(State is Loyal) ,', beta)), at = .75,
     tick = F, pos = .52)


x.points3 <- lowess(x$alpha[x$group == 3], x$beta[x$group == 3], f = .05)$x
y.points3 <- lowess(x$alpha[x$group == 3], x$beta[x$group == 3], f = .05)$y

x.points2 <- lowess(x$alpha[x$group == 2], x$beta[x$group == 2], f = .05)$x
y.points2 <- lowess(x$alpha[x$group == 2], x$beta[x$group == 2], f = .05)$y

polygon(c(x.points2, rev(x.points2)), c(y.points2, rep(.5, length(y.points2))), col = 'lightgrey')
polygon(c(x.points3, rev(x.points3)), c(y.points3, rep(.5, length(y.points3))), col = 'darkgrey')

lines(lowess(x$alpha[x$group == 2], x$beta[x$group == 2], f = .05))
lines(lowess(x$alpha[x$group == 3], x$beta[x$group == 3], f = .05))
abline(0,1, lty = 3)

points(.8, .7, pch = 18, col = 'black')

## plot legend


plot(0,
     col = 'white',
     bty = 'n',
     xaxt = 'n', yaxt = 'n')
legend('topleft',
       title = expression(paste('Removal Threshold, ', hat(b), ', in:')),
       y.intersp = 1,
       pch = c(NA, NA, NA, NA, NA, 18, NA),
       lty = c(NA, NA, NA, NA, 3, NA, NA),
       lwd = c(NA, NA, NA, NA, 1, NA, NA),
       fill = c('white', 'lightgrey', 'darkgrey', 'white', 'white', 'white', 'white'),
       border = c('black', 'black', 'black', 'white', 'white', 'white', 'white', 'white'),
       legend = c('Self-enforcing region',
                  'Enforceable region', 'Unenforceable region', '', '45-degree line',
                  expression(paste(alpha, ' and ', beta, ' used in Figure 1' )),  '',
                  expression(paste(sigma, '=1, ',
                                   phi, '=1, ', 'k=1'))), cex = 1)




dev.off()

























